%  Copyright (C) 2002 Regents of the University of Michigan, portions used with permission 
%  For more information, see http://csem.engin.umich.edu/tools/swmf
\chapter{Equations and Physics \label{chapter:equations}}

\section{The MHD Equations \label{section:mhd}}
The \BATSRUS\ code solves the equations of
magnetohydrodynamics, or MHD.  The current distribution of \BATSRUS\
uses the ideal MHD equations, meaning that the pressure is a scalar,
there is no heat flow and the conductivity is taken to be infinite
(magnetic resistivity is zero).  In this limit the MHD equations can be
written as follows.  We give both the so called primitive and
conservative forms.

The primitive forms is:
\begin{equation}
\label{eq:continuity_primative}
\frac{\partial \rho}{\partial t} + \nabla\cdot \left(\rho {\bf u}\right)
= 0
\end{equation}
\begin{equation}
\label{eq:momentum_primative}
\rho \frac{\partial {\bf u}}{\partial t} + \rho\left({\bf u}\cdot
\nabla\right){\bf u} + \nabla \cdot \left[\left( p + \frac{B^2}{2
\mu_0}\right)I
-\frac{\bf BB}{\mu_0}\right] = \rho {\bf g}
\end{equation}
\begin{equation}
\label{eq:bfield_primative}
\frac{\partial {\bf B}}{\partial t} + \nabla\cdot\left({\bf u\,B} - {\bf
B\,u}\right) = 0
\end{equation}
\begin{equation}
\label{eq:pressure_primative}
\frac{1}{\gamma-1} \frac{\partial p}{\partial t} + \frac{1}{\gamma-1}
\left({\bf u}\cdot \nabla\right)p + \frac{\gamma}{\gamma-1}
p \left(\nabla \cdot {\bf u}\right) = 0
\end{equation}

The conservative form is:
\begin{equation}
\label{eq:continuity_conserve}
\frac{\partial \rho}{\partial t} + \nabla\cdot \left(\rho {\bf u}\right)
= 0
\end{equation}
\begin{equation}
\label{eq:momentum_conserve}
\frac{\partial \left(\rho {\bf u}\right)}{\partial t} +
\nabla\cdot\left[
\rho{\bf u\,u}+
\left(p + \frac{B^2}{2 \mu_0}\right)I
- \frac{{\bf B\,B}}{\mu_0}\right] =
\rho {\bf g}
\end{equation}
\begin{equation}
\label{eq:bfield_conserve}
\frac{\partial {\bf B}}{\partial t} + \nabla\cdot\left({\bf u\,B} - {\bf
B\,u}\right) = 0
\end{equation}
\begin{equation}
\label{eq:energy_conserve}
\frac{\partial\varepsilon}{\partial t}
+  \nabla \cdot \left[ {\bf u} \left( \varepsilon + p
+ \frac{B^2}{2 \mu_0} \right) - \frac{\left({\bf u}\cdot{\bf B}\right) {\bf
B}}{\mu_0}\right] =
\rho {\bf g}\cdot {\bf u}
\end{equation}
where $\rho$ is the plasma mass density, $\bf u$ the plasma velocity, 
$\bf B$ the magnetic field, $p$ the pressure, ${\bf g}$ the
gravitational acceleration due to the central body and  $\varepsilon$
is the total energy density given by 
\begin{equation}
\varepsilon=\frac{\rho u^2}{2}+\frac{p}{\gamma-1}+\frac{B^2}{2\mu_0}
\end{equation}


\section{Normalization \label{section:normalization}}

In the \BATSRUS\ code, the above MHD equations are solved in their
normalized form.  This is done so that constants like $\mu_0$
are not lost in the coding of algorithms.  It also helps
to ensure that numbers are near order 1.  For the
MHD equations, we need only three independent factors for normalizing.
We have chosen to use a reference density ($\rho_0$), 
a reference velocity ($u_0$) and the
typical length scale of the problem ($R_0$).  The normalizations of
the typical quantities are: 
\begin{alignat}{4}
\tilde{t}     &= \frac{u_0}{R_0} t                         & \qquad \qquad
\tilde{\bf r} &= \frac{\bf r}{R_0}                         & \qquad \qquad
\tilde{\rho}  &= \frac{\rho}{\rho_0}                       & \qquad \qquad
\tilde{\bf u} &= \frac{\bf u}{u_0}                           \nonumber \\
\tilde{p}     &= \frac{p}{\rho_0 u_0^2}                    & 
\tilde{m}     &= \frac{m}{\rho_0 R_0^3}                    &
\tilde{\bf B} &= \frac{\bf B}{\sqrt{\mu_0 \rho_0 u_0^2}}   &  
\tilde{\bf E} &= \frac{\bf E}{u_0^2 \sqrt{\mu_0 \rho_0}}   \nonumber \\
\tilde{\bf j} &= \frac{R_0 \sqrt{\mu_0}}{\sqrt{ \rho_0 u_0^2}}{\bf j} 
\end{alignat}
where the variables with a tilde (\~\ ) are unitless, or normalized.

In general, the upstream solar wind values for density ($\rho_0 =
\rho_{SW}$), sound speed ($u_0=a_SW$), and planet radius ($R_0 = R_P$)
are chosen as the normalizing parameters. 


\section{Gravity \label{section:gravity}}

Gravity is included and computed in the code assuming a
spherical body.  Gravitational forces are computed for each cell, but
include only forces due to the central body.  This means
that there is no ``self gravitation'' of the plasma. Gravity can
be turned off for each individual case by changing the way that 
the initial conditions are loaded in {\tt setICs.f90}.  Note that
for Earth gravity is currently turned off.

\section{Mass Loading \label{section:mass_loading}}

Mass loading is the process of adding new plasma to the system.  This
mass addition usually results from ionization processes or through
charge exchange.  Since these are non-MHD processes, they cannot be
dealt with directly.  We can, however, include their effects
by adding source terms to the above MHD equations.  Since these terms
are problem dependent, we will let the following suffice as an example.
If we are ionizing a neutral gas which is at rest, the mass addition
rate would be $\nu \rho_n$, where $\nu$ is the ionization rate and
$\rho_n$ the neutral mass density.  The MHD equations, in primitive form,
for this process are the following
\begin{equation}
\label{eq:continuity_primative_mass_loading}
\frac{\partial \rho}{\partial t} + \nabla\cdot \left(\rho {\bf u}\right)
= \nu \rho_n
\end{equation}
\begin{equation}
\label{eq:momentum_primative_mass_loading}
\rho \frac{\partial {\bf u}}{\partial t} + \rho\left({\bf u}\cdot
\nabla\right){\bf u} + \nabla \cdot \left[\left( p + \frac{B^2}{2
\mu_0}\right)I
-\frac{\bf BB}{\mu_0}\right]
 = -\nu \rho_n {\bf u}
\end{equation}
\begin{equation}
\label{eq:pressure_primative_mass_loading}
\frac{1}{\gamma-1} \frac{\partial p}{\partial t} + \frac{1}{\gamma-1}
\left({\bf u}\cdot \nabla\right)p + \frac{\gamma}{\gamma-1}
p \left(\nabla \cdot {\bf u}\right) = \frac{1}{2}\nu \rho_n u^2
\end{equation}
where the equation for magnetic field has not been included because it
is unchanged.  The terms on the right hand side are the source terms 
and can be included in the code to model the ionization and pickup
process.  The equations above are given in the primitive form, because for the
case that we have given, the momentum and energy density 
source terms are zero when the equations
are converted to the conservative form.

\section{Other Source Terms \label{section:other_sources}}

Mass loading is not the only physical process that can be modeled
through source terms.  \BATSRUS\ is designed to use any source term
that the user desires.  Heating and divergence $\bf B$ control (see
section~\ref{section:divb}) are examples of other source terms.  Since the source terms are
not used in determining the time step, there are issues of stability
when the source term dominates the solution in a region.

\section{Intrinsic Magnetic Field \label{section:B0}}

Because many solar system bodies posses intrinsic magnetic fields,
\BATSRUS\ allows the user to specify such a field.  Currently, 
a multipole expansion of the field from dipole up to octupole terms 
is available, but extension to higher multipoles would be trivial.  The
magnetic dipole axis can be tilted and can move as a function
of time.  The dipole tilt occurs only in the x-z plane of the
simulation where the x-axis points at the sun and the z-axis is chosen
to contain the magnetic field dipole axis.  This means that we are in a
Geocentric Solar Magnetic (GSM) coordinate system.
Future implementations will allow the user to specify a
more complicated magnetic field model such as the International
Geophysical Reference Model (IGRF) for Earth.

Typical planetary intrinsic fields have steep gradients which can
present problems for numerical models.  In order to alleviate this
problem, \BATSRUS\ splits the intrinsic magnetic field from the 
rest of the magnetic field.
\begin{equation}
{\bf B} = {\bf B}_0 + {\bf B}_1
\end{equation}
Here ${\bf B}$ is the total magnetic field, ${\bf B}_0$ is the
intrinsic magnetic field and  ${\bf B}_1$ is the deviation of the field
from ${\bf B}_)$.  Note that ${\bf B}_1$ is not required to be small and
that when derivations are done using the split field no terms are
dropped.  
The only restrictions are that $\nabla \cdot {\bf B}_0 = 0$ and
$\nabla \times {\bf B}_0 = 0$. 
The split is analytic, but allows \BATSRUS\ to advance only
${\bf B}_1$ with time and thereby avoid the problems of the steep
gradients in ${\bf B}_0$.

As noted above, the total magnetic field is split into the intrinsic part 
and the deviative part.  Since we allow $\bf B_0$ to change with time, 
there should be a $\frac{\partial {\bf B}_0}{\partial t}$ term in the 
induction equation.  Terms proportional to this would also appear in
other equations.  This term is currently left out, but will be included
in a future update.


\section{Rotation \label{section:rotation}}
Rotation of the central body is implemented through boundary conditions.
As with the magnetic field, the user can specify the rotation axis.
Currently the rotation axis can only be in the x-z plane.  This means
that one component of the rotation axis (the one perpendicular to the
x-z plane) is currently ignored.  This will be changed 
in future distributions.

\section{Boris Correction \label{section:boris}}

For the standard MHD equations, equations~(\ref{eq:continuity_primative})--
(\ref{eq:energy_conserve}), the Alfv\'en speed is
defined as
\begin{equation}
v_A = \frac{\bf B}{\sqrt{\mu_0 \rho}}
\end{equation}
This is not limited to the speed of light and can in fact exceed the 
speed of light in locations where the magnetic field is strong (near the
poles of the  Earth, Jupiter and Saturn, for example).  This fact has large
implications in determining the time step which the code uses.  Since the 
time step in the numerical simulation is determined by the fastest wave speed, 
it is the Alfv\'en speed which often limits the size of a time step.

In order to resolve this physically inadmissible speed, Boris pointed out 
that if one retains
the displacement current term in Amp\`eres law, and then derives the
MHD equations, a solution where the Alfv\'en speed is
limited to the speed of light is obtained.  While this solution is completely 
physical, it still allows velocities that are near the speed of light,
which in turn can still produce a very small time step.  

With the Alfv\'en speed limited to the speed of light, we can now 
use a method to increase the time step in the code.  If one were
to artificially lower the speed of light, the Alfv\'en speed would
also be limited to this new lower value and the time step in the
code would be increased.
Making this change to the speed of light is a common technique in 
the numerical MHD community.

Both the Boris correction (inclusion of the displacement current term)
and the ability to artificially lower the speed of light are implemented 
in \BATSRUS.  We should point out to the user that the term ``Boris
correction'', when used by the MHD community, usually means both the 
limiting of the Alfv\'en speed to the speed of light by including the
displacement current term and the artificial lowering of the speed of
light.  While the
including of the displacement current and the 
corresponding limiting of the Alfv\'en speed to the speed of light is 
completely physical, the lowering of the speed of light is not.
At this point we should warn the user that while lowering the speed of
light is a widely used
technique, it has not been well investigated.  We know that limiting
the speed of light limits the wave speeds and therefore we believe that this
will change the solution.  The option to lower the speed of light should be used
with care.


\input{ionosphere.tex}






